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ROBUST STABILITY OF LINEAR SYSTEMS - SOME 
COMPUTATIONAL CONSIDERATIONS* 

by 

Alan J. Laub** 


1, INTRODUCTION 

In this paper we shall concentrate on some of the computational issues 
which arise in studying the robust stability of linear systems. Insofar 
as possible, we shall use notation consistent with Stein's paper [1] and 
we shall make frequent reference to that work. 

As we saw in [11 a basic stability question for a linear time -invariant 
system with transfer matrix G(s) is the following; given that a nominal 
closed- loop feedback system is stable, does the feedback system remain 
stable when subjected to perturbations and how large can those perturba- 
tions be? It turned out, through invocation of the Nyqxiist Criterion, 
that the size of the allowable perturbations was related to the "nearness 
to singularity" of the return difference matrix I + G(jci)} , Closed-loop 
stability was said to be "robust" if G could tolerate considerable 
perturbation before I + G became singular. 
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research was partially supported by NASA under grant NGL-22-009-124 and 
the Department of Energy under grant ET-7S- (01-3395) . 

Laboratory for Information and Decision Systems, Rm. 35-331, M. I. T., 
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We shall now indulge in a modicum of abstraction and attempt to 
formalize the notion of robustness. The definition will employ some 
jargon from algebraic geometry and will be applicable to a variety of 
situations. While no deep results from algebraic geometry need be em- 
ployed, the exercise of formulating a precise definition is a useful one 

for clarifying one's thinking. 

N 

Let p € IR be a vector of parameters from some problem being studied 

and suppose we are interested in some property H of this data. The vector 

p may consist of the elements of various matrices, for example. If IT 

is true at some nominal parameter set p^ we are frequently concerned with 

whether H' remains true in a "neighborhood" of p^. 

Por example, p_ may be the elements (a,,, ..., a, , a.,,..., a } 

0 11 In 21 nn 

of a nonsingular nxn matrix and we are interested in the nonsingularity 
of nearby matrices . We shall proceed to formalize the often-heard statement 
that "almost all nxn matrices are nonsingular". First, the jargon; 

f N T 

Definition 1 ; A variety 1/ = {p 6 IR ; • • • ' ^ “ If , k) 

where (x^, . . . , x^) e 3i[Xj^, . . . , x^^l are polynomials . 

V is proper if t/ ^ and nontrivial if V 

Definition 2 ; A property is a function IT; IR^^ {o, l} . The property 
n holds if HCp) = 1 and- fails if II (p) * Q, 

Definition 3 t If 1/ is a proper variety, H is generic relative to 1/ 
provided H (p) =0 only if p S 1/. A property H is 
generic if such a 1/ exists. 

Our discussion to this point is purely algebraic. Mow let us intro- 
duce a topology on say the topology induced by some vector norm [| • [[ . 
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rtxrthennore, let 1/ be any nontrivial, proper variety. Then we have 
the following topological definition. 

Definition 4 ; The property II is well-posed at p 6 1/ (the complement of 
V) if H also holds in a sufficiently small neighborhood 
of p. 

Lemma 1 ; The set S of points where a generic property is well-posed 

Q 

is open and dense. Moreover, the Lebesgue measure of S 
is zero. 

The proof of Lemma 1 is routine and is omitted. It is easy to see 
that a point p where a generic property holds is well-posed but that the 
converse is not necessarily true. 

We now have sufficient framework to make a formal definition of 
robustness. 

Definition St Given a point p with generic property II (generic with 

respect to some proper variety [/) well-posed at p, let 

d “ min [{ d - v [[ . 
ve(/ 

We say H is‘ robust at p if d is "large". 

The number d is frequently difficult to compute or estimate, tflhen 
it can be determined, it gives valuable information about how much 
perturbation or uncertainty can be tolerated at p. For the situation 
of special interest in this paper, Example 2 below, we shall see that 
d can be explicitly calculated, at least theoretically. We now illustrate 
the above concepts with two examples. 
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Example 1 

This exas5»le is chosen from Wonham [21 who uses the concepts of 
genericity and well-posedness in nontrivial ways for a variety of control- 
theoretic problems. In this trivial example, we seek solutions of the 
system of linear equations 

Ax B b 

where A (i.e., A is an nixn matrix with real coefficients) and b SIR™. 

Our parameter vector is p where 

T N 

p — 1 ^ * f t » * • t a f b_,.««, b ) € IR f N — mn m 

i-L in mn 1 m 

T 

( denotes transpose) , II is the property of the equation having a solution 

which is equivalent/ of course, to the statements that b S Im A or 

(1 2 \ 

rk[A, b] = rk A. For example, if A and b 

(0 if b, 2b, 

11(1,2,2,4; b ,b,) = ) ^ i 

• I 1 if b^ = 2bj^ 

It is then easy to show the following: (see [2]) 

1. H is generic if and only if m < n. 

2, H is well-posed at p if and only if rk A = m. 

Example 2 

This example is similar to Example 1 in the special case m - n. We 

are given a nonsingular matrix A € and we are concerned with the 

T 

nearness of A to singularity. Identi^ing A with p = (a^^^/ . . ., 

a ) we define the property II by 
2x nn 
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0 if p represents a singular matrix 


II(P) » < 


if p represents a nonsingular matrix . 

Then it is easy to see that II is a generic property and well-posed where 
it holds. This Is the precise statement that "almost all nxn matrices 
are nonsingular". Formally writing down the determinant of A as a poly- 
nomial in a, a defines the necessary variety I/. It turns out, 
lx nn 

« 

in a theorem attributed by Kahan [3] to Gastinel, that the distance d 

c 

from a point pel/ to 1/ can be explicitly determined. 


Theorem 1 ; A nonsingular matrix A differs from a singular matrix by no 
more in norm them — i- — i.e., given A, 

rNi 


— ■ min{|lE II ! A + E is singular} 

IIa II 


Thus d = = and we might say that A is robust with respect to 

|a" I 


invertibility if d is "large". To avoid certain scaling difficulties, 
it may be more desirable to work with a relative measure of distance, 
d^®^, defined by 


^rel _ d 1 

11a 11 IIaII • 11a"^1 


1 

K(A) 


The quantity K(A) is recognizable as the condition number of A with 
respect to inversion. Of course, all the above quantities depend on the 
particular matrix norm used. To exhibit the specific dependence on the 
norm |1 • we shall append a subscript "q". For example. 





d - — i — . 

The ninimizing E in Theorem 1 can be explicitly constructed for a nttmber 
of standard matrix norms. For example: 

Let A have singular value decomposition A » USv'^ where U, V 6 

are orthogonal and S » diag{a , a }. The a. 's, a, a >0, 

J. n 1 1 — ’ — n 

m 

are the singulau: values of A. The minimizing E is given by E = URV 
where R» diag{0,..,, 0, -c }. Then 



1 


and A + E is singular. The singular direction, i.e., a nonzero 
vector z such that (A + E)z = 0, is given by the n— column of 


n 

Halloo” Jt'sx { J] I a , . j } , {l,2,..., n} . 

i6n j=l 

-1 -1 ” 

Suppose A « [ct . ] and )| A |j = I I k 6 n. Then the 

minimizing E is a matrix all of whose elements are 0 except 
tH 

for the k-— column which consists of the elements 


- =9" “ki “to 



In fact, letting z = sgn 


with the only 
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3. 


th T 

nonzero component of u being in the k— - row, we have E » -zu 


and clearly [js|| 


^ . Now, (1+SA ^)z » (1 - u*a”^z)z“ 0 


since the k— element of is f l“k-}l “ II^'*^iloo®° u’^a'^z - 1. 

j-1 


-1 


Hence A+E«(I + EA )Ais singular. Moreover, the singular direction is 


given by a“^z since (A+E)a“^z » 0. 


n 


A [L - max (I |a I>. 
^ jen i-1 


The results for this norm are analogous to li * i| „ and can be derived 
directly or by noticing that || A j| j^ ■ [| . For completeness we 


note that if |I A |L ■ k 6 n and 

J- IK 



then the minimizing E is given by E 



We shall see in Section 3 how the results in Example 2 can be applied 
in studying robustness of stability of linear systems. 
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2. THE LINEAR SYSTEMS SETTING 

In this section we shall provide a brief introduction to both the 
linear tine-invariant systenis setting 2 uid to the fundanental notion 
of feedback. This will serve a two-fold purpose: first, to set the stage 
for the basic steU^llity results ^uld second, to introduce the jargon and 
notation, especially for non-engineers. This material is standard and 
can be fotxnd in any of a number of standard textbooks on control systems. 

We shall consider modelling physical systems by models which take 
the form of a system of linear constant-coefficient ordinary differential 
equations 


x(t) 

= Ax(t) + Bu(t) 

(1) 

y(t) 

» Cx(t) 

(2) 


The vector x is an n-vector of states, u is an m-vector of inputs or 
controls, and y is an r-vector of outputs or observed varieibles. 

Starting from the initial condition x(0) the solution of (1) is 
well-known to be 


x(t) = 


e x(0} 




(t- T)A 


‘Bu(T)dT, t > 0 


(3) 


SO that the output is given by 


y(t) = Ce^\(0) + he 


oh 


(t - T) A, 


Bu(T)dT 


tA . 


where e is the matrix exponential defined, 
computed, by 


tA r t A 

' ~ V o 

k ==0 


, t ^ 0 

but not generally 


( 4 ) 
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The matrix is called the i mpulse response matrix. 

Denoting (one-sided) Laplace transforms by upper case letters, taj<e 
Laplace transforms in (4) to get 

y(s) - CX(s) - C(SI- A)“^x(0) + C(sl- A)”^BU(s) . (5) 

The matrix G(s): ■ C(sl - A) is called the tr^msfer matrix. Notice 
that G(s) is the Laplace transform of the impulse response matrix. 

As will M seen in the sequel, it is of interest to study the 
response of the above linear system to sinusoidal inputs of the form 

u(t) « t ^ 0 (6) 

where v is a constant m-vector, (O is the frequency of the sinusoidal 
input, and j = »/^. The response of (1) to this input can then be 
shown to be of the form 

x(t) « e^\ + (jO)I- A)"^Bve^^^, t^O (7) 

where a is a constant n-vector depending on initial conditions. Now, 

in the case where A is stable (i.e., its spectrum lies in the left- 

tA 

half of the complex plane) the quantity e a goes to zero as t approaches 
+«. The resulting output 

y(t) = C(jo)I - A) ^Bve^^^ (B) 

is called the steady-state frequency response and the matrix 

G(jco) ; = C(jwl - A)"S , (9) 

which turns out to be the transfer function evaluated at s = jw, is 
called the frequency response matrix. 
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Turning now to the case of a real signal given by 


Uj^(t) » Vj^sin(ajt + tjj) » t ^0 (10) 

u^(t) ■ 0/ i =» 1, . , . , m; i k, 

we have steady-state frequency response of the Zth output given by 


y^(t) » Ivj^ sin(0Jt + (11) 

where » arg(G^j^(jaj)) . 

Aside from its obvious importance in the above analysis, the 
frequency response matrix is important for two reasons: 

1. Sinusoidal signals are readily available as test signals 

for a linear system so G(jw) can be experimentally determined. 


2. Various plots or graphs associated with G(j6j) can be used to 

analyze control systems, for exan^le, with respect to stability. 
Plots such as those associated with the names of Bode, Nichols, 
and Nyquist are essentially different ways of graphically 
representing arg(G^j^(jo)) ) as functions of 

< 0 . These plots are used extensively in the analysis of 
single-input single-output control systems where the robust- 
ness of stability, e.g., the amount of gain and phase margin 
available, is checked essentially visually. The appropriate 
techniques in the multiple-input multiple -output case are 
still being investigated and part of the motivation for the 
research in [1] and this paper is directed towards this end. 
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Turning now to the notion of feedback whose essential Idea Is to 
allow for stability of a system in the face of uncertainty (noise, 
model error, etc.), the diagram below illustrates the basic (unity) 
feedback control system: 



Pig. 1. Basic Feedback Control System 

Here u is a reference input, y is the output, and e =» u - y is the error 
or difference between the reference input and the output which we wish 
to be, ideally, zero. The plant, compensators, actuators, and sensors 
are all represented by G. There are much more elaborate and detailed 
feedback structures than that described above and the structure can be 
studied in a considerably more general function-space setting (see (4) , 
for example) than the simple linear causal time-invariant setting we 
shall consider. However, the simple system is adequate to exhibit most 
of the key ideas in this paper. Now, in this system we have 

e = u- y = u-Ge (12) 

or. 


(I+G)e = u 


(13) 




- 12 - 


The quantity I + G is called the return dijeference matrix. As in [1] , 

the matrix G(jai) then provides sufficient data, via the Nyquist criterion, 

to test for stability of the closed-loop system. Henceforth, we shall 

assume that our nominal feedback system above is stable in which case 

* 

I+G is invertible. Then from (13) we have 

e - (1 + G)"^u (14) 

so that 

y ■ Ge « G(I + G) . (15) 

In (IS), the quantity G(s) (I + G(s)) ^ is called the closed-loop transfer 
matrix while G(jco) (I + G(j(o)) ^ is called the closed-loop frequency 
response matrix. We then pose the basic stability question; 

Does the nominal feedback system remain stable subjected 

to perturbations and how large can those perturbations be? 

Let us observe at this point that there is nothing sacred about 
linearity in the above discussion and more general nonlinear treat- 
ments can be found in [4] and [S] , for example. The question of "near- 
ness to singularity" of (1 + G), even in the nonlinear case, is naturally 
intimately related to a notion of condition number for nonlinear 
equations. The interested reader could readily adapt the ideas of 
Pheinboldt [6] to the particular application at hand here. 
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3. BASIC STABILITY RESULTS AND RELATED TOPICS 
a. A DDITIVE AND MULTIPLICATIVE PERTURBATIONS 
We shall consider two fundamental types of perturbations in the 
basic feedback system of Pig. 1. Throughout this section, || • [j will 
denote any matrix norm with |[ 1 fi “ 1. The first case to be considered 
is tb.e case of additive perturbations to G, pictured belowj 



Pig. 2. Additive Perturbations 

In other words, the nominal G is perturbed to G -f L. Under the assumptions 
that both the nominal closed-loop system and the perturbation L are 
stable it can be seen from the Nyquist criterion and the identity 

I + G + L = (I + G) [I + (I + G)'^L] (16) 

that the perturbed closed-loop system remains stable if 

[j(I + G(jOJ))“^L(j6)) [[ < 1, CO > 0 (17) 

A weaker condition than (17) but one which directly exposes L is 






-14- 


||L(jOJ) i| < - — i — , 0) > 0 (18) 

|I(I + G{j(u)) -^11 

The second case to be considered is that of multiplicative perturba- 
tions : 



Fig, 3. Multiplicative Perturbations 

In this case, the nominal g is perturbed to G(I+L) . Under the assumptions 
that both the nominal closed-loop system and the perturbation L are 
stable it can be shown from the Nyquist criterion and the identity 

I + G(I+L) = (I + G) [1+ (I + g"^)“^L] (19) 

that the perturbed closed-loop system remains stable if 

|[(I + G’^(j£D))'^L(ja3) li <1, 6) > 0 (20) 

(assuming G ^ exists). Again, a weaker condition than (20) but one 
which directly exposes L is 

liL(ja»|l < i ^,W>0 . (21) 

|(l + G"-"(jW)}'-"[ 
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Remark 1 ; As we noted in Section the above inequalities are tight, 
i 4 e , , the < cannot be replaced with ^ . 

Remark 2; Where convenient we shall henceforth drop the "jw" arguments 
Remark 3; It must be stressed that the results based on 

11(1+ II 1|l II < 1 (18) , (21) 

'are weaker than those based on 

|i(l + G~^)"^L II < 1 (17), -(20) 

since 

k 

^ ||(i+ g“S"^l|| 1 ild + G'^-r^II • 1 |l||, (22) 

% 

+1 II 

For example, if L = cd + G" ) for some constant c, |c| < 1, the 
differences in the bounds are obvious. In (13), (21) we have 

|l(l+'G"S“^|i • IIlII = |c|*K(I +G“^) 
while in (17), (20) we have 



and it is possible to have 

Icj « |c| • K(l+G“^) . 

However, for random perturbations L, (22) is often approximately an 
equality. To see this, note that a random (dense) L will almost surely 
be invertible; recall Example 2. It is then easy to show that 
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Again, since L is random« it will aln^ost surely be well-conditioned 
(w.r.t. inversion) so that ||l ^|| * • Hence, 

A related aspect, also worth noting, follows from the inequalities 

< ||(I + G^S-l|| • IlLlI . 

K<I+ G"-^) 

+1 

If (1 + G ) is reasonably well-conditioned (K(I + G ) near 1) , the 
majorization (22) will not be a bad overestimate. 

Remark 4; By our discussion in section 1, the appropriate measure of 
stability robustness is 

d - min . ( 23 ) 

03>0 !|(I + G-^(j(0))“'^|| 

and in the sequel we shall consider methods of efficiently plotting 

^ - — z — as a function of o). This quantity is familiar from 

classical sensitivity analysis where it is shown, in the single-input 
single-output case, that the change in the output of a closed-loop 
system, due to (additive) perturbations in G (scalar) , is reduced by 
a factor of 1 + G compared with the open-loop effect. 

Remark 5; So far we have required nothing of our norm other than 
II I II = 1. Of course, a frequently occurring norm in much of the 
analysis of linear systems is the spectral norm jl'H, . 


In that case 
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1 ■‘‘1 

■ j ' Si — smallest singular value of (I+G~ ). Let 

— Xfc "A M 


11(1 + G=^) -•‘I 


dg(m) 


li + s 


-"(3<«))‘"ll, 


( 24 ) 


We are Interested in plotting d (<o) versus (o for large numbers of 

<I 

ca's. We shall see in the sequel that determining <^2 

somewhat more expensive to determine than, say d, (io) or d (w) . More- 

X ® 

over, note that 


^lUIL < (25) 

vm 

and 

^ lUl! < ||a|[^ < ^ 11a [Ij (26) 

vm 

for A € Since we are usually most interested in order-of-magnitude 

estimates of d (oi) , d_ (w) will lie in a strip sufficiently close to 
dj^(co), for example, to give the same qualitative information. The 
number m which is the number of inputs/outputs in the system is typically 
no more than about 10 and is frequently much less. 

b. RELATIONSHIPS BETI-?EEN ADDITIVE AND MULTIPLICATIVE PERTURBATIONS 

The following theorem relates additive and multiplicative perturba- 
tions. Again, the "jCd's" will be omitted for convenience and all 
relations will be assumed to hold for all co > 0 . 
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Theoren 2: 


|{(1 + g“^)"^|| - i|(I+G)‘-^|| 


- 1 | 


< 1 


Proof; From the identity 


(I + g"^)"^ + (1+ G)"^ = I 


(27) 


we have 

||[(I+g“^)“^[| - [l(I + G)'^|(|£|j(I+G“^)"^ + (1 + G)"^|I = III II » 1, 
We now get immediately the following useful corollary *• 

Corollary It Assuming that both the nominal closed-loop feedback 
system of Pig, 1 and the perturbation L are stable then the perturbed 
system is stable under; 


(a) additive perturbations if 

1 


L < 


1 + (i + g‘^)"^ 


(28) 


(b) multiplicative perturbations if 


L < 


1 + (I + G) 


-li 


(29) 


Proof ! Follows immediately from Theorem 2 noting that 


1 + j|(I + G“-^) 




c. SPECIAL RESULTS FOR THE SPECTRAL NORM 

In this subsection we shall present some results related to those 
in subsections a. and b. but specisdized to the ][• - norm. For 
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a matrix H e C*™*® with singular values cr, (H) > ...> a (H) > 0 we note 

X ~ “ m 


-li 


that II H j}^ * If H is nonsingular, ||H Ij^ 

II • II2 “ norm (28) becomes 




> 0. In the 


aj^(L) < 


a (i + g“^) 

m 

1 + a (1 + g“^) 

m 


while (29) becomes 


a (i + G) 


m 


1 + cr (1 + G) 
n 


We shall make great use in the sequel of the following result 
of Pan [7] . 


Theorem 3: Let A,B S Then 

(a) o,. . (A+B) < a. (A) + a. (B); i > 1, j > 1 

l+3“X — 1 3 — — 

< a. (A)a. (B); i > 1, j > 1 

Part (b) of Theorem 3 can be used to relate 0 (I + G) and 0 (1+G ^). 

m m 

Theorem 4 ! (a) — ^ — 0 (I+G*^) <0 (I + G) < || G ||.o (I + G*^) 

Ib'-Il2 " - " - ' '2 m 

(b) — i a (I + G) <0 (i + G*^) < IIg*^ IL a (I + G) 

IIgIIj ” - - ‘ "2 " 

Proof; Follows immediately from Theorem 3 using 

1 + g"^ = g’^(I + G) and I + G = G(I + g”^) . 
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For the rest of this subsection we shall let H denote either 

-1 

I + G or I + G according to whether additive or multiplicative 
perturbations are appropriate. The next theorem will show how the 
singular values of H -f-. L can be bounded in terms of j| L and the 
singular values of H. 


Theorem 5; Suppose ^ > 0 for some k, 1 ^ 5. W/ and 

IIlIIj < S. Suppose further that 8 < Then; 


(a) 




h"^l) 


1 - J- 

“k 


(b) + L) > - 0 . 

(Note; If k m, H + L is not necessarily invertible if S is too large.) 


•L •*!. ^2* '*•1 

Proof; (a) Use ISI + K L-H L and A » I + H L, B =» -H L, i « k, 

j a m-k + 1 in Theorem 3(a) to get 

a (I) < a (I + h”^l) + a . , (h“^l) . 
ro ■— Jc m— ic+i 

-1 -1 

Thus a (I + H L) > 1 - a , (H L) 
k ~ m-k+1 

>_ 1 - II L II 2 ' (h"^) by Theorem 3 (b) 

= 1 “ 1(l ll2'^3,(H) 


(b) Use H F H + L - L and A = H + L, B = -L, i = k, j = 1 
in Theorem 3 (a) to get 

L) > a^CH) - ||l H 2 > - 8 . 
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The case k » nj is of special interest in Theorem 5 as it bears 
directly on our two basic inequalities (18) and (21) of the form 



which are sufficient to guarantee stability of a perturbed closed-loop 
system. Specifically, if || h”^|| j 1 ^ II ^ II 1 3 with 0 j< 0 < a, 
then H+L is invertible and |1(H+L)"^|| £ or a^(H+L) - 0. 

Note that Theorem 5 was expressed in terms of isolating || L jj2 . 3y 
analogy with the inequalities (17) and (20) we can also have the fol- 
lowing stronger, but perhaps less useful, theorem. 

Theorem 6 : Suppose 1 - <5 where 0 < 5 < 1 and 1 < k < m. 

Then; 

(a) + h“^L) ^ 5 

(b) a. (H + L) > — 

Proof ; - (a) From the proof of Theorem 5 we have 
a. (i + h"^l) > 1 - a , (h"^l) > 6 

k — m-iC+l — 

(b) From I + H = H ^(H+ L) and Theorem 3(b) we have 
+ h"^l) + L) . IIh’^II^ 

whence a, (H+L) > — — ^ . 

^ “ IIh'^IIs 
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d. SPECIAL RESULTS WHEN G(s) » C(sl - aTS 


In this subsection we shall nuUce use of the fact 
response matrix is of the form 

G(jai) - C(jwi - A)"^B 

Let us further define 

that the frequency 

F(jW) ■ C(j6Jl - A + BC)"^B 
Recall the Sherman -Morrison-Woodbury formula: 

(W + XYZ)*^ = W’^ - w"^X(y"^ + 2W"^X)“^ZW“^ 

(30) 

assuming the indicated inverses exist. Then it is easy to verify that 

-1 

(I + G(jai)) ■ = 1 - p(jOJ) 
and# from (27) , 

-1 -1 

(31) 

(I + G (jOJ)) = F(jw) 

(32) 

Thus our results in the last section (for example 
and 6) can all be cast in terms of p by noting that 

, Theorems 4, 5, 

and 

(33) 

Moreover, 

(34) 

||(I + G)-"|| = 

(35) 
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and 

i|(l+G-l)-’-|| - ||f|I (36) 

for any of the norm we have been considering (in particular, k " m in 
(33) and (34)). Use of (31) and (32) results in an apparent savings 
in the number of linear systems to be solved (i.e., number of inversions) 
emd we shall exploit this fact in the next section. 
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4. COMPUTATIONAL PROBLEMS 

a. COMPUTATION OF FREQUSNCy RESPONSE MATRICES 

As we have seen above, an object of considerable interest in studying 

the robustness of staOsility of linear systems is a graph of - ■ " ' ■ . 

11(1 + G“^(j(o))‘^!| 

as a function of OJ. When Gtjw) ■ C{jWI - A) we saw that [|(1 + G(j6j))~^|j » 
[Il-P(jtu)|l and |((I + G"^(jco))”^{| » |{p(jco)(| where P{j(u) » C(jwi - A + BC)"^B, 
Thus, regardless of the norm used, a q\iantity of the form 

C(j(0I - K)“^B (37) 

must first be computed. We shall assume throughout this and the next sub- 
section that; (i> B C H are given 

(ii) n > m 

(iii) (37) is to be evaluated for a large number, N, of 
values of W; typically N » n. 

Rather than concentrate on exact operation counts, which may be fairly 
meaningless anyway, we shall give only order-of-magnitude estimates. 

It will be seen that the bulk of the computational load rests on evalu- 
ating matrices of the form (37) and so we shall focus initially on 
that problem. 

nxn _1 

If A 6 3R is dense, the most efficient evaluation of C(jtoI - A) B 

by an LU factorization of A, solution of m triangular systems to get 
-1 

(jwi - A) B, and finally a matrix multiplication, requires approximately 
13 12 2 

yn + ymn + m n multiplications (and a like number of additions; we 
shall henceforth count only multiplications). This figure, when multiplied 
by N, represents a rather large amount of computation. 
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If A Is initially transformed, however, the computational burden 
ceui be reduced quite considerably. If T is a similarity transformation 
on A we have 

C(jwl - A)“S = CT{j«I - t"^AT)“^t"S . 

Let us define 

H ■ t”^at 

^md agree, for convenience to still label CT, T the transformed C and 
B matrices, respectively, as C, B respectively. We now have the problem 
of evaluating 

C(j«I - H)“^B 

where H may now be in such a form that (jcoi - H) ^ can be computed in 

3 

less than 0(n ) operations. For example, A can always be reduced to 

5 3 

upper Hessenberg form by (stabilized) elementary transformations {— n 

5 3 

multiplications) or by orthogonal transformations (y n multiplications) . 

3 

These trsuis formations are very stable numerically and, while 0(n ), are 

performed only once at the beginning of the calculations. The resulting 

linear system to be solved - for N different values of (O — now has an 

upper Hessenberg coefficient matrix and can be solved in approximately 
1 2 

— mn multiplications. Moreover, Hessenberg systems czm be solved very 

accurately with the growth factor in Gaussian elimination bounded above 

by n; see [8] . Computing C(jO)l - H) ^B still requires an additional 
2 

m n multiplications. Neglecting the initial transformation and deter- 
mination of CT and T ^B, the Hessenberg method requires approximately 
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12 2 

~ rtm + n n multiplications (for each value of co), a considerable savings 

3 

over the 0(n ) algorithm if n » m. 

Of course^ other trems formations T are possible. One possibility is 
to reduce A to upper triangular (Schur) form by means of orthogonal simi- 
larities. This is considerably more expensive than reduction to upper 
Hessenberg but, again, need only be done once at the beginning. How- 
ever, the resulting linear system to be solved at each step is upper 

2 

triangular and so still requires 0(mn } multiplications. Because of 

potential difficulties with multiple eigenvalues of A there seems to be 

little real advantage gained by this procedure. Substantial savings 

could be gained though if the eigenstructure* of A were such that it . 

was diagonalizable by a reliably computable T. Since this involves 

consideration of the essentially open numerical problems associated with 

computing invariant subspaces we shall not pursue the details here. 

But assuming such a transformation were possible, C(j(i3l - D) with 

2 

D diagonal, could be computed with approximately mn + m n multiplications 

for each value of to. Attractive as this appears, the potential for severe 

ill-conditioning of the eigenproblem associated with A render this latter 

method unreliable as a general-purpose approach. We shall subsequently 

consider only the Hessenberg method. 

The analysis above has been done under the assumption that complex 

-1 

arithmetic was performed. We now outline how G = C(jcoi - H) B might 
be determined using only real arithmetic. The matrix H is assumed to 
be in upper Hessenberg form. We wish to solve first 


(jcol - H)Z = 3 


(38) 
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Then 

G ** C2 • 

Suppose Z « X + jY where X, Y ejFP”". Upon equating real and imaginary 
paurts in (38) we get the following order 2n real system to determine 
X and Y: 



Thus X » i HY and Y - -w(£o^l + The matrix (co^I + H“) will be 

2 2 

invertible if (jcoi - K) is invertible. Note that (w I + H ) is no longer 
upper Hessenberg but is almost in the sense of having two rather than one 
nonzero subdiagonal. Its shape is wholly typified for n = 5 by the 
matrix 



Linear systems involving matrices of this type can be solved using 
2 

approximately n multiplications. We summarize the Hessenberg method 
using real arithmetic: 


(i) Reduce A to upper Hessenberg form H, transform B and C, 

2 

and conpute H ; this step is done only once. 
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2 2 

(ii) Solve (Co I + H )Y » -coB for Y, 

(iii) Compute X ■ ~ HY . 

(iv) Compute G » (CX) + j (CY) . 

2 

step (ii) requires approximately mn multiplications, step (iii) requires 

12 2 
approximately — mn , and step (iii) approximately ai sj. The total number 

3 2 2 

of multiplications is approximately ^ mn + m n . 

Storage requirements for the Hessenberg method with real arithmetic 
are approximately double those for complex arithmetic. 

b. COMPUTATION OF ROBUSTNESS MEASURES 

We have seen above that quantities of the form (37) can be reliably 
2 

evaluated in 0(mn ) operations. There then remains the problem of 
determining (35) or (36). 

Case 1 ; [j • II2 

For (35) , the singular value decomposition (SVD) of I + G(j<o) can 

be computed for each value of co. Each SVD typically requires approximately 
3 

6m multiplications. The smallest singular value is then the quantity of 
interest. For (36) , inversion of G can be avoided by finding the SVD 
of F(jCij), again in approximately 6m^ multiplications. The inverse of 
the largest singular value of ? is then the quantity of interest. 

Case 2 ; i| • ||^ or || • [| ^ 

Use of either of these norms in (35) or (36) involves negligible 

2 

computation as compared to Case 1, namely about m additions and absolute 
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values and tn'-l arithmetic comparisons. 

In both cases, the additional work required is usually small compared 
2 

with 0(mn ) especially if n » m. However, if m is large relative to 
n, significant savings can be realized in using ||* ||^or rather than 
[j * j[ 2 > In fact, using our previous approximate operation counts for 
the Hessenberg method amd setting n =* ten,' we have 

work per value of Oi using |j • | [j ~ k^ + 2k + 12 

P a I 

work per value of 6) using j| ■ |[^ or || • |(^ k + 2k 

2 

k + 2k + 24 

Note though that p ~ - if singulcir directions are also com- 

k + 2k 

puted. 

In the event A (or A - BC) can be successfully diagonalized as 
mentioned in Section 4. a. the potential savings in avoiding j{ • are 
somewhat greater. In fact, we then have 

k+6 

k+12 

(or p St — ^ — if singular directions are also computed) . 

The above coit 5 >arisons are only approximate and should in no way 
be construed as definitive statements. The purpose of this section is 
to merely introduce certain aspects of the numerical computations and 
suggest further avenues of exploration. A great deal of numerical 
experimentation remains to be done. Reliable software such as LINPACK 
[9] for linear systems will be of great benefit in this research. 
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5. CONCLUSIONS 

We began this paper with an attempt at a "formal" definition of 

robustness. We then applied the definition to the problem of robustness 

of stability of linear systems as discussed in [1] . The cases of both 

additive and multiplicative perturbations were discussed and a number of 

relationships between the two cases were given. Finally, a number of 

computational aspects of the theory were discussed including a proposed 

new method for evaluating general transfer or frequency response matrices 

The new method is numerically stable and efficient, requiring only 
2 

0(mn ) operations to update for new values of the frequency parameter 
rather than O(n^) 

A number of interesting research areas suggest themselves in this 
work. One such area is that of constradned perturbations. For example, 
in our basic problem we were concerned with the nearness to singularity 
of a nonsingular matrix A 6 If the admissible perturbations E 

are somehow constrained for one reason or another, for example E upper 
tricingular, the usual bound on || E || for which A + E is singular but E 
is "dense" may be overly pessimistic. Related to this is the fact that 
our bounds were derived for the "worst case". The size of perturbations 
allowed in a linear system to ensure continued closed- loop stability may 
very well be larger than we have derived if inputs to the system are 
constrained in certain directions. 

We have concentrated in this paper on the analysis of linear control 
systems . There are many interesting — and difficult — synthesis problems 
however. For example, can A, B, C be chosen to assign certain singular 

4*1 

values of I + G ? What is the effect of changes in B or C on the 
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behavior of I.I + G ? Can a matrix K be determined so that I + (GK}~ 
has certain singular values? 

On the computational side, more research needs to be done on updating 
parametric problems. That is, suppose we have a matrix (say, G(joj)) 
which depends "in a rank m way" on a parameter w. When (o changes how 
can various quantities be updated efficiently? 

Finally, as mentioned in Section 4.b., a great deal of numerical 
experimentation is necessary to get a qualitative feel for the numbers 
in determining robustness measures. 
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